home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / ACOSH.C next >
Encoding:
C/C++ Source or Header  |  1989-03-21  |  3.0 KB  |  156 lines
  1. /*                            acosh.c
  2.  *
  3.  *    Inverse hyperbolic cosine
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, acosh();
  10.  *
  11.  * y = acosh( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns inverse hyperbolic cosine of argument.
  18.  *
  19.  * If 1 <= x < 1.5, a rational approximation
  20.  *
  21.  *    sqrt(z) * P(z)/Q(z)
  22.  *
  23.  * where z = x-1, is used.  Otherwise,
  24.  *
  25.  * acosh(x)  =  log( x + sqrt( (x-1)(x+1) ).
  26.  *
  27.  *
  28.  *
  29.  * ACCURACY:
  30.  *
  31.  *                      Relative error:
  32.  * arithmetic   domain     # trials      peak         rms
  33.  *    DEC       1,3         10000       4.1e-17     1.1e-17
  34.  *    IEEE      1,3         30000       4.6e-16     8.7e-17
  35.  *
  36.  *
  37.  * ERROR MESSAGES:
  38.  *
  39.  *   message         condition      value returned
  40.  * acosh domain       |x| < 1            0.0
  41.  *
  42.  */
  43.  
  44. /*                            acosh.c    */
  45.  
  46. /*
  47. Cephes Math Library Release 2.1:  December, 1988
  48. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  49. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50. */
  51.  
  52.  
  53. /* acosh(z) = sqrt(x) * R(x), z = x + 1, interval 0 < x < 0.5 */
  54.  
  55. #include "mconf.h"
  56.  
  57. #ifdef UNK
  58. static double P[] = {
  59.  1.18801130533544501356E2,
  60.  3.94726656571334401102E3,
  61.  3.43989375926195455866E4,
  62.  1.08102874834699867335E5,
  63.  1.10855947270161294369E5
  64. };
  65. static double Q[] = {
  66. /* 1.00000000000000000000E0,*/
  67.  1.86145380837903397292E2,
  68.  4.15352677227719831579E3,
  69.  2.97683430363289370382E4,
  70.  8.29725251988426222434E4,
  71.  7.83869920495893927727E4
  72. };
  73. #endif
  74.  
  75. #ifdef DEC
  76. static short P[] = {
  77. 0041755,0115055,0144002,0146444,
  78. 0043166,0132103,0155150,0150302,
  79. 0044006,0057360,0003021,0162753,
  80. 0044323,0021557,0175225,0056253,
  81. 0044330,0101771,0040046,0006636
  82. };
  83. static short Q[] = {
  84. /*0040200,0000000,0000000,0000000,*/
  85. 0042072,0022467,0126670,0041232,
  86. 0043201,0146066,0152142,0034015,
  87. 0043750,0110257,0121165,0026100,
  88. 0044242,0007103,0034667,0033173,
  89. 0044231,0014576,0175573,0017472
  90. };
  91. #endif
  92.  
  93. #ifdef IBMPC
  94. static short P[] = {
  95. 0x59a4,0xb900,0xb345,0x405d,
  96. 0x1a18,0x7b4d,0xd688,0x40ae,
  97. 0x3cbd,0x00c2,0xcbde,0x40e0,
  98. 0xab95,0xff52,0x646d,0x40fa,
  99. 0xc1b4,0x2804,0x107f,0x40fb
  100. };
  101. static short Q[] = {
  102. /*0x0000,0x0000,0x0000,0x3ff0,*/
  103. 0x0853,0xf5b7,0x44a6,0x4067,
  104. 0x4702,0xda8c,0x3986,0x40b0,
  105. 0xa588,0xf44e,0x1215,0x40dd,
  106. 0xe6cf,0x6736,0x41c8,0x40f4,
  107. 0x63e7,0xdf6f,0x232f,0x40f3
  108. };
  109. #endif
  110.  
  111. #ifdef MIEEE
  112. static short P[] = {
  113. 0x405d,0xb345,0xb900,0x59a4,
  114. 0x40ae,0xd688,0x7b4d,0x1a18,
  115. 0x40e0,0xcbde,0x00c2,0x3cbd,
  116. 0x40fa,0x646d,0xff52,0xab95,
  117. 0x40fb,0x107f,0x2804,0xc1b4
  118. };
  119. static short Q[] = {
  120. 0x4067,0x44a6,0xf5b7,0x0853,
  121. 0x40b0,0x3986,0xda8c,0x4702,
  122. 0x40dd,0x1215,0xf44e,0xa588,
  123. 0x40f4,0x41c8,0x6736,0xe6cf,
  124. 0x40f3,0x232f,0xdf6f,0x63e7,
  125. };
  126. #endif
  127.  
  128. extern double LOGE2;
  129.  
  130. double acosh(x)
  131. double x;
  132. {
  133. double a, z;
  134. double log(), sqrt(), polevl(), p1evl();
  135.  
  136. if( x < 1.0 )
  137.     {
  138.     mtherr( "acosh", DOMAIN );
  139.     return(0.0);
  140.     }
  141.  
  142. if( x > 1.0e8 )
  143.     return( log(x) + LOGE2 );
  144.  
  145. z = x - 1.0;
  146.  
  147. if( z < 0.5 )
  148.     {
  149.     a = sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5) );
  150.     return( a );
  151.     }
  152.  
  153. a = sqrt( z*(x+1.0) );
  154. return( log(x + a) );
  155. }
  156.